dir_data <-"../Data/lidar"
year <- "2013"
dirname <- paste(dir_data,"vector",year,"treetops",sep="/")
ws <- 18
test_filename <- paste0("treetops_lmf_ws",ws)
test_filename_full <- paste(dirname,test_filename,sep = "/")
# filename <- "../Data/lidar/vector/2013/treetops/treetops_lmf_ws18"

dir_crowns <- paste(dir_data,"raster","crowns",sep="/")
dir.create(dir_crowns,recursive = T)

## Warning in dir.create(dir_crowns, recursive = T): '../Data/lidar/raster/crowns'
## already exists

## Read underlying raster
raster_filename <- paste(dir_data,"raster","raster2013.tif",sep="/")
r <- raster(raster_filename)

treetops <- readOGR(test_filename_full)

## OGR data source with driver: GeoJSON 
## Source: "/home/lune/Data/ai4er/mres/lidar/vector/2013/treetops/treetops_lmf_ws18", layer: "2013_ws18"
## with 2510 features
## It has 2 fields

plot(r)
plot(treetops,add=T)

filename_dalponte <- paste(dir_crowns,
                           paste0(
                             paste(year,ws,"dalponte",
                                sep = "_"),
                             ".tif"
                           ),
                           sep="/")
# crowns_dalponte <- dalponte2016(
#   r,
#   treetops,
#   th_tree = 2,
#   th_seed = 0.45,
#   th_cr = 0.55,
#   max_cr = 50,
#   ID = "treeID"
# )()
# writeRaster(crowns_dalponte,filename_dalponte)
crowns_dalponte <- raster(filename_dalponte)

plot(crowns_dalponte,col=pastel.colors(1000))

filename_silva <- paste(dir_crowns,
                           paste0(
                             paste(year,ws,"silva",
                                sep = "_"),
                             ".tif"
                           ),
                           sep="/")
# crowns_silva <- silva2016(r,treetops)()
# writeRaster(crowns_silva,filename_silva)
crowns_silva <- raster(filename_silva)

plot(crowns_silva,col=pastel.colors(1000))

# contours <- rasterToPolygons(crowns_dalponte,dissolve = TRUE)

subsampling <- 100000
pixel_area <- 0.5*0.5
ncells <- r@ncols*r@nrows


hist <- hist(crowns_dalponte,breaks=2510,plot=F, maxpixels=subsampling)

## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 1%
## of the raster cells were used. 100000 values used.

areas <- sort(hist$counts)*ncells/subsampling*pixel_area
plot(areas,pch=".")

hist <- hist(crowns_silva,breaks=2510,plot=F, maxpixels=subsampling)

## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 1%
## of the raster cells were used. 100000 values used.

areas <- sort(hist$counts)*ncells/subsampling*pixel_area
points(areas,pch=".",add=T,col="red")

## Warning in plot.xy(xy.coords(x, y), type = type, ...): "add" is not a graphical
## parameter